i#-*-encoding:latin-1-*-
# Disponivel no pacote de programas como: slice.py
from numpy import *
from scipy.stats import *

class Amostrador:
    def __init__(self, alvo, proposta):
        """
        Incializa o amostrador
        alvo: dados definindo densidade alvo
        proposta: objeto gerador de propostas.
        """
        self.alvo = kde.gaussian_kde(alvo)
        self.prop = proposta
    
    def Run(self,n,burnin=100):
        '''
        Roda o amostrador por n passos
        '''
        amostra = zeros(n,float)
        alvo = self.alvo
        i = 0
        while i < n-1:
            x = amostra[i]
            y = uniform.rvs(0,alvo.evaluate(x))
            inova = self.prop.rvs()[0]
            can = x + inova
            aprob = min([1,alvo.evaluate(can)/alvo.evaluate(x)])
            u = uniform.rvs()[0]
            if u < aprob:
                amostra[i+1] = x = can
                i += 1
        return amostra[burnin:]
        
if __name__=="__main__":
    import pylab as P
    dados = concatenate((norm.rvs(size=500),norm.rvs(4,1,size=500)))
    a = Amostrador(dados,norm(0,1))
    res = a.Run(10000)
    P.plot(arange(-5,10,.01),a.alvo.evaluate(arange(-5,10,.01)),'r',lw=2)
    P.hist(res,normed=1,alpha=0.5)
    P.legend(['Alvo','Amostra'])
    P.show()
